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Abstract 

The aim of this paper is to investigate a reaction-diffusion Leslie-Gower predator-prey 
model, incorporating the intraguild predation and both self and cross-diffusion. The 
longtime behaviour of the solutions is analysed, proving the existence of an absorbing 
set. The existence of patterns is investigated by looking for conditions guaranteeing 
that an equilibrium, stable in the absence of diffusion, becomes unstable when diffusion 
is allowed. 
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1 Introduction 


The analysis of the occurrence of diffusion-driven instability fits into the theoretical 
framework of reaction-diffusion models. The use of these models exhibits a particular 
efficacy in fields ranging from morphogenetic processes [1, 2], to oscillating chemical 
reactions [3, 4], and population dynamics including epidemiological models [5, 6] 
and predator-prey models [7, 8]. The model analysed in this work is of predator-prey 
type. According to the pioneering works of Lotka and Volterra [9, 10], the classical 
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predator-prey models describe the interaction between two populations in which a 
species (the predators) feeds on the other (the preys). In order to better describe some 
real situations, various generalizations to the classical model have been proposed. An 
example is given by food chain systems which are a multi-level trophic model [11] 
and predation is a relevant kind of relationship between the levels. Food chain models 
often describe three trophic-level webs, composed of logistic prey and Lotka- Volterra 
or Holling type II specialist predator and top-predator [12]. The model we study here is 
especially based on a modified version of the Leslie-Gower scheme [13-15], that has 
been extensively studied in the last decades [16, 17, and references therein]. A Leslie— 
Gower model is a predator-prey model where the carrying capacity of the predator 
is proportional to the number of prey [15]. Modified versions of the Leslie-Gower 
predator—prey model with various functional responses have been deeply investigated 
(see [18, 19, and references therein]), in particular, in [18] the authors analyzed a 
Leslie—Gower predator-prey model in which the top-predator population feeds only 
on the meso-predator population and influences its consumption rate, focusing the 
analysis on how the top-predator interference in the capture of the meso-predator 
affects the populations’ dynamics. 

Usually, the predator—prey interactions are modelled by functional responses, which 
depict how the consumption rate of individual consumers changes with respect to 
resource density and describe ideal populations interactions neglecting spatial hetero- 
geneity. According to Fick’s law, the species are spatially heterogeneous and hence 
individuals will tend to migrate toward the regions of the lower population to add the 
possibility of survival [20]. The aim of this paper is to generalize the model proposed 
in [18], supposing that the populations are heterogeneously mixed in the environment. 
To do that, we added to the system linear self and cross-diffusion terms, in order to 
analyse the occurrence of Turing instability. 

The paper is organized as follows. In Sect.2 the mathematical model is presented. 
In Sect. 3 we proved the boundedness of the solutions of the system and the existence 
of an absorbing set in the phase space. In Sect. 4 the biologically meaningful equilibria 
are determined and discussed. Sections 5 and 6 deal with the linear instability of the 
coexistence equilibrium, with the aim to determine conditions which guarantee the 
stability in the absence of diffusion and diffusion-driven instability, namely Turing 
instability. Section7 discusses the numerical investigations on the stability results we 
found. The paper ends with a concluding Section that recaps all the obtained results. 


2 Mathematical model 


Let us consider three interacting species in a bounded domain Q C RY (N = 2,3), 
whose evolution is described by a predator-prey model. Let us denote by P(X, f), 
M(X, í) and T(X, 7) the prey, meso-predator and top-predator population density, 
respectively. The model describes a tritrophic web where a prey population serves 
as the only food for the meso-predators. This population, in turn, serves as a favourite 
food for the top-predator, hence both populations meso and top predator are spe- 
cialists. Following [18], we employ a Tanner model [21] to describe the interaction 
dynamics between top-predator and meso-predator and a Gause model [22] to describe 
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the interference of top-predator on the interaction dynamics between prey and meso- 
predator. According to [23], let us suppose that if the population is homogeneously 
mixed in the environment: (i) the prey population grows logistically in absence of the 
meso-predator population with an intrinsic growth rate 7 and environmental carry- 
ing capacity K; (ii) the functional response of the meso predator to the resource is 
of Beddington-De Angelis type and is given by ar where the term cT in the 
denominators describes the interference of top-predators on meso-predators; (iii) the 
functional response of the top- Pedara, describing the meso-top predator dynamics, is 
an Holling-Type II given by 777 Ti M ; (iv) the top-predator growth equation is of logistic- 
type with a modification of the conventional one. The environmental carrying capacity 
of the top-predator population is not constant but proportional to meso-predator abun- 
dance, i.e. it is equal to mM, where the parameter m is the measure of the food quality 
that the meso-predator provides for conversion into meso-predators births. To gener- 
alize the results obtained in [18], we assume that the populations can spread randomly 
in the environment. This assumption is useful for modelling, for example, the typical 
behaviour of predators that move in search of prey or of prey that escape from high-risk 
predation zones. Therefore, the mathematical model is 


oP S P aM D 
—~ = r{l1 ~ a +yAP 
ot K 1+bP+4+cT 


aM aqP eT ~ 5 
— =M = - - d| + PRAM (1) 
ot 1+bP+¢T 1+7M 
ai ST|1 3 + RAT 
— =5 2 
ot mM 3 
where y;, (i = 1, 2,3) are the constant diffusion coefficients, while all the other 


biological pare: are assumed to be non-negative. In the sequel, we will denote 


by A = {K,a,b,¢, d, č, ñ, ñ, q, F, 5, Y1, Y2, 3} € RM. To system (1) we append 
the ere Neumann (no-flux) boundary ndon 


VP-n=0, VM-n=0, VT-n=0, ondQx Rt (2) 


with 0Q the boundary of Q and n the outward unit normal of 0&2, and the following 
smooth positive initial conditions 


P(x, 0) = Po(x), M(x,0) = Mo(x), T(K,0)=To(X), XER. (3) 
Moreover, let us suppose that: (i) Q is of class C? (p > 2), with the interior cone 


property; (ii) p € W!?(Q) N W!7(dQ), Vi > 0, Vy € {P, M, T}. 
Introducing the following dimensionless parameters 


a a a e 
=l _ å — aq E — db — ¢mr emb 
E da 8 Da Oe P= a 8 ae O 
y yı y2 y ¥3 
l= 3793 A = Soa 9 3 = 275> 

FL? 7L2 FL? 
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L being the diameter of Q, the initial boundary value problem (1)—(3) becomes 


ðu 1 v PIA 
— =u u u, 
ot b+u+cw vı 


ðv u ew FERE (5) 
aaa v 
ot q” b+u+cw n+uv raD 


ðw (1 x) BA 
— = ws { 1 — — w, 
ot v É 


Vu-n=0, Vv-n=0, Vw-n=0, onðQ x Rt (6) 
u(x, 0) = uo(x), v(x, 0) = v(x), wx, 0) = wo(x), xe Q. (7) 
In the sequel, we denote by (-,-) and ||-|| the scalar product and norm in L?(Q), 


respectively, by ||-||oo the norm in L% (Q) and by A= {b,c,d,e,n,g,5, Vi} € R10 
the set of the dimensionless biological parameters. 


3 Longtime behavior of solutions 


This Section deals with the boundedness of the solutions of system (5)—(7). Let T>0 
be an arbitrary fixed time and let us denote by Rg = Nx (0, 7] the parabolic cylinder. 
Defining T} = dQ x [0, 7), the parabolic boundary of Q; is Fe =T,U{Qx {t= 
O}}, the following theorems hold. 


Theorem 3.1 Let (u, v, w) € [C? (27) N COD be anon negative solution of (5)— 
(7). Then Yọ € {u, v, w}, ọ is bounded almost everywhere in Q according to 


u(x, t) < max fi, max wc :=M1, v(x, t) <Coo(vo(X)) := M2, 
Q 


(8) 
w(x, t) < max {Mo max w0] := M3 
Q 


where Coo is a positive constant depending on the initial data. 


Proof Let us set max u = u(X1, t1). We distinguish the following cases. 
T 


1. If (x1, t1) € Qz, then (5); implies that 


ou 
| yı^u—u(l J <0 (9) 
at (x14) 


Since SH Nt) = 0, Au|(x,,4,) < 0, (9) holds if u(x), t1) < 1. 
2. If (x1,t1) € Q x {t = 0} then u(x;, 0) < max uo(x). 
Q 
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3. If (x1, 4) € Ds, then by virtue of the regularity of the domain &, since Q verifies 
in any point y € ðX the interior ball condition, there exists an open ball B* C Q 
with y € 0B*. If x; € dQ and u(x, tı) > 1, choosing the radius of B* sufficiently 

a 


small, it follows that [yı Au — $F ly, > 0 for each y € B*. By virtue of Hopf’s 


Lemma [24], one obtains that du l&i) > 0, which is not possible in view of the 
boundary conditions (6). 


Hence u(x, t) < max{1, max ug(x)}. In order to prove (8)2, in view of (5) it turns out 


that u(x, t) is a sub-solution of the problem 


ox 
VX -n=0 ondQ~x Rt (10) 


X(x, t) = Xo(x) = max vo(x). 
Q 


Since qX < x + i; in view of Theorem 1 of [25], one obtains that, denoting by 
Tt(Xo) the maximal existence time of the solution X(x, t) of (10), since - from the 
continuous dependence on the initial data - there exists a positive constant C(X0) 
such that || X(-, t)|| < C(Xo0), Vt € (0, tT(Xo)), the solution X (x, t) exists for all time 
and there exists a positive constant Coo such that || X(-, t)lloo < Coo(X0), Vt > 0, 
and (8)2 is proved. 

In view of (8)2, following, step by step, the same procedure used to prove (8), 
inequality (8)3 can be obtained. o 


Theorem 3.2 Ve > 0 the manifold 
X; = fon v, w) € (Rt)? : lull? + full? + lwl? < a roe (11) 
with 
ā=2qd, b=2I9| [a +dq)M? + qM2 + (s + qd) M3] (12) 


in an absorbing set for system (5). 


Proof Let us set E(t) = llul? + lvl? + ||w||?. Multiplying (5) by (u, v, w)’, inte- 
grating over Q and adding the resulting equations, by virtue of the divergence theorem 
and the boundary conditions (6), it turns out that 


l dE 2 2 2 

5 ap SAGE + + ga) lull? + alol? + (s + ga) wl. 
In view of (8), one obtains that dE < —aE + b. Following step by step the procedure 
in [26], one can prove that X, is an absorbing set. o 
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4 Biologically meaningful equilibria: existence and a priori estimates 


The biologically meaningful equilibria are the non-negative solutions of the system 


v 
A 1 = 0, 
” | 7 | 
u ew 
A d| =0, (13) 
y2 vta] n+v | 


yz Aw + ws (1- À =0. 


Theorem 4.1 Let (u(x), v(x), w(x)) € [C7(Q) N C!(Q)P be a positive solution of 
(13) subject to zero flux boundary conditions. Then, setting 0 < y < min{y1, 2, 3}, 
there exist three positive constants Cj(y), (i = 1, 2, 3) depending on A and Q such 
that: 


maxu(x) <1, maxv(x) < ay = M2, maxw(x) < Mp, 
a u(x) i maxa v(x) maxa w(x) Cs) 
— 2 < Cy), —2— < Cly), —2— < 67) 
ming u(x) ming v(x) ming w(x) 
Proof Let us set š} € Q such that u(X1) = maxu(x). By virtue of the max- 


imum principle (see [24]), one gets —Au(x,) > 0. Hence, from (13); one has 


l-u- > 0, which implies (14),. Regarding inequality (14), let us 


v 
b+u+cw A 


consider S(x) = yıqu(x) + y2u(x) and X2 € Q such that S(x2) = max S(x). Since 
Q 


~AS(X) > 0, one obtains [ua =p ee dv]. > 0, that, by virtue of (14), 
%2 


implies v(X2) < 1/d, i.e. 
1 
min v(x) < —. (15) 
Q d 


Inequalities (14)4—(14)6 are recovered by using the Harnack’s inequality (see (27, 28]). 

From (15) and (14)s, (14)2 is immediately obtained. Finally, setting x3 € Q such that 

w(x3) = max w(x) and applying the maximum principle, from (13)3 and (14)2 one 
Q 


obtains w(X3) < v(%3) < M2. Oo 
Ag =-— in Q 

Let us set a; the lowest positive eigenvalue of the spectral problem £ aa 
Vo-n=0 ondQ 


and @ = a Ís gdQ Vo € {u, v, w}. The following theorem provides sufficient con- 
ditions for the non-existence of non—constant solutions of (13). 
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Theorem 4.2 If 


Cy) ð 
> 1 
nZ a (i + oR ) 


Ca(y)M2 [(qcūù eq(in+d) eqw sdū 
> 16 
a ai 2b? 2n2 n? + 20 (16) 
X C3(y)M2 Acv F cqu 4 eq(n + v) pdw s e abtaqcw , 
ay 2b? 2b? 2n? 20 M b+ut+cw 


holds, then the system (13) does not admit any positive non-constant solution. 


Proof Let (u, v, w) be a positive solution of (13). Multiplying (13); by A(u — ū)/u, 
(13)2 by (v — v)/v and (13)3 by (w — w)/w, integrating over Q and adding the 
resulting equations, by virtue of the divergence theorem and boundary conditions (6), 
one obtains 


tots i. 3 Ù 5 
Ayı | Vu d+ y | Vv d+ ys | — Vwa 
Qu Q Uv Q Ww 


= (—bA— Aū — Act) + b i eee? A DA 
= HSRC taD ET Je ú 
-g2 Se 
+a f Ca = d2 acò | La l 
Qg (b +u +cw)(b +u + cw) o (b +u +cw)(b +u + cw) 


af (v — v)(w — w) d2- ( 4 Df @ — V)(w — w) jio 
cqu en v E S E E 
1 |. bEukem beipen) ET EET 


7 (v—v) sw f v—d)\(w— wD) (w — Ñ)? 
+eqiv yt z L i dQ sf m . a7 


qb + qcù 
b+itcw’ 


By virtue of the Poincaré inequality and of Theorem (3.1), choosing A = 


from (17) one gets 


a [ u — Dd + < f on oe I (w — w)2dQ 
AGM C2(y) M2 C3(y) M2 
qcūü eq(n+v) eqw sw 
<A dQ 
Gaze = 2b? i) Í (uid ( 2b? ü 2n? j n? T =) 
e Acūù cqū eq(n+0) Sw s Í - 
2 2 
—v)dQ = dQ 
x fe ee (ia Tr Dn Doum =) oe 
(18) 
that is impossible when (16) holds. oO 


For the sake of completeness, in the sequel the results of [18] regarding the existence 
of meaningful equilibria are recalled. 


Theorem 4.3 A necessary condition for the existence of a generic equilibrium to (13) 
isu <l. 
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Theorem 4.4 Ifd > 1, there are no equilibrium points. 


Avoiding the null solution, the boundary equilibrium is the top-predator free equilib- 


rium E; = (24, p048 0), existing only ifd < = [18]. 


The interior equilibrium Æ = (i, 0, W) with non-null components represents the 

coexistence of all three populations. Let us remark that, since {u > 0, v > 0, w > O}, 

(13); admits a positive root only if u < 1. From (13)z, one has that: 
ü 


e since = =< > 
b+ut+cw 1+b 
u 


there exists a positive solution only if d < 


1 

1+b’ 

u i nA ‘ EN bd 
- — < —, there exists a positive solution only if u > —. 
(b+u+cw) b+u l-d 


Then, the condition 


e since 


d < —— (19) 


is necessary for the existence of the coexistence equilibrium and, if E exists then, one 


bd n 
necessarily has T <u < l, v= w. The coexistence equilibrium F exists if 


n 1 : u<1l (20) 
max )7_ 4" Ps (eas 
and is given by 
b+u)—u 
jets es (21) 
l+cu-—c 


being u a positive root of the equation 
Aut + Bu? + Cu? + Du + E = 0, (22) 
where 


A=c, B=1-2c+be—c? —d-e, 

C=-1+b-—c—2be + 2c? +d —2bd +e —2be+cdn, 

ee eee eee) E ee ee eee ee a 
E=b'd+b’e+bdn — bedn. 


Let us remark that 
c>1, b(e+d) <dn(c—-1), (24) 
are sufficient conditions such that (22) admits at least one positive root u. 
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5 Preliminaries to stability 


Let us denote by (u*, v*, w*) a generic constant equilibrium. Introducing the pertur- 
bation fields {U =u — u*, V = v — v*, W = w — v*}, system (5) becomes 


aU ~ 
=ayU+ay2V +ay3W+yAU + Fi 
aV ~ 
gp TU + a22V +anW + AV + Fo (25) 
aw i 
aE = a31U + a32 V +433W + AW + F 
where 
1 2ut v*(b+cw*) | u* ; cu*v* 
nes “ (b + u* + cw*)2’ S b+ u* + cw*’ m pau Ee] 
qv*(b + cw*) — qu* enqw* f 
ars (b + u* + cw*)2’ maa= b+u*+cw* (n+v*)? dq; (26) 
oF * *\2 * 
423 = a. ai a31 =0; a32 a an=s{ a 
v 


; =s : 
b+u*+cw* n+v* (vx)? 


and the non-linear terms F(U, V, W), (U, V, W), R (U, V, W) are defined as 
follows: 


z ela (QV + v*)(b + c(63W + w*) 2 2c? (01U + u*)(@2V + v*) 2 
I= (b +01U +u* + c(03W + w*))3 (b+ OU +u* +c(03W + w*))3 
b+c0 * 8 j 
c(@3W + w*) UV 4 c(0,U + u*) VW 
(b+ 6,U + u* + c(03W + w*))2 (b+6,U + u* + c(03W + w*))2 
( c(02V + v*) ORLAR AAR v 
(b+6,U +u* +c(03W +w*))?  (b+01U +u* + c(03W + w*))3 
= 2q (02V + v*)(b + c(03W + w*)) 5 2enq(03W + w*) 2 
(b +01U +u* + c(03W + w*))3 (n+ 0V + v*)3 (27) 
2q OU + u*)(62V + v*) j 
| (b+ 0U +u* +c(03W + w*))3 
+ * * 
q(b+ c(63W + w*) UV cq(0,U + u*) VW 
(b+ 0,U +u* + c(63W + w*))? (b+ 0,U +u* + c(03W + w*))2 
1 (b+ c(0 * 
cq(nV +1”) lg TE w 
(b+6,U +u* +c(03W +w*))2?  (b+0U +u* + c(03W + w*))3 
AK 2s(03W + w*)? 5 2s 2, 2s(03W + w*) 
(02V + v*)3 V +0* (02V + v*)2 


with 6; € (0,1), (i = 1,2,3). To system (25), the following no-flux boundary 
conditions are associated 


VU-n=0, VV-n=0, VW-n=0, onðQ xR. (28) 
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6 Linear instability analysis 


In this Section, we will investigate the linear instability of the coexistence equilibrium 
E. In particular, we will look both for conditions guaranteeing the stability in absence 
of diffusion and for conditions guaranteeing the diffusion-driven instability, namely 
Turing instability. Let us recall that in [18] the linear stability conditions for the 
boundary equilibrium Æ; in absence of diffusion have been determined and it has 
been proved that the boundary equilibrium is linearly unstable. 


6.1 Linear instability in absence of diffusion 


In absence of diffusion, i.e. for yi = 0 (i = 1, 2, 3), the linear system associated to 
system (25) evaluated around E is 


ax 411 412 413 
a £°X, with X=(U,V,W), L} = | an an az |, (29) 
0 s =s 


where X is the state vector and £? is the Jacobian matrix evaluated in the coexistence 
equilibrium Æ, according to (26). The characteristic equation associated to £° is 43 — 
I ba + Ev — I = 0, where IP (i = 1, 2, 3) are the principal invariants in absence of 
diffusion. The interior equilibrium E is linearly stable if and only if the Routh-Hurwitz 
conditions are verified [29-31]: 


P20, R20; PRP-R<O0. (30) 


Conditions (30) obviously imply IR > 0. If at least one of (30) is reversed, there 
exists at least one eigenvalue of £? with positive real parts, and hence, in that case, 
the linear instability of E is guaranteed. We remark that, in view of (26) evaluated in 
E , conditions 


ayt< 0, a22 < 0, (31) 


are sufficient to guarantee the validity of (30). In order to derive more biologically 
significant conditions, let us look for simpler sufficient conditions guaranteeing (30). 
Setting 


the following theorem holds. 
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Theorem 6.1 If 


b(e+d) 1 
d< —, c>l, n>—W—, m>-e, 
1+b d(c — 1) 2, (33) 
d —d)[n+ v — m) _ (1 — d)[n + v2(1 — m)? 
e> , m<u<il ; 
nvı(l — m) env] 


then the coexistence equilibrium E is linearly stable. 


Proof We recall that condition (33); is necessary for the existence of E, while (33)2- 

(33)3 are sufficient for its existence. Furthermore, if E exists, then (20) holds true. 
5(b ja 

From (31), it follows that aj; = 1—2u T mn < 1 — 2u. By virtue of (20), 

condition (33)4 guarantees that aj; < 0. In view of (31) and (20), it turns out that 


— ee . From (20) and (21) and by virtue of (32), it easily turns 
(n+ 0)? 
envi (l — a) | 


[n + v2(1 — m)]? 
Condition (33) implies that a22 < 0, while (33)5 guarantees that u > m. The thesis 
is obtained since it has been proved that (33) implies a;; < 0, Vi = 1, 2. oO 


an <q|l-—d 


out that vı (1 — u) < v < v2(1 — u). Hence az < q [i 


6.2 Self-diffusion driven linear instability 


ax 
The linear system associated to (25) is a = LX + DAX, where 


a11 412 413 7190 0 
X=(U,V,W), £= |an an az], D=|0%0 (34) 
0 s =s 007 


The dispersion relation governing the eigenvalues À in terms of wavenumber k is [30] 
P= TEKNA? + Bn(k?)A — hô) = 0, (35) 


where 


Tk (k?) = (LÀ) — KD) = 19 —P n+ +73), 


hi?) =Khiintuntny)—PIn(an—s)+2(a11-s) +73 (a11 +a22)]4 R, 


2 6 4 (36) 
h(k?) = —k?det(D) + k [a11 y2V3 + a2aviy3 — $V1 V2] 
Kiyas — azs) — vrais + y3 (411422 — 412421)] + 1, 
All roots of (35) have negative real part if and only if 
Ty(k*) <0, h(k?) <0, Ty(k?)In(k*) — h(k?) <0, Vk?. (37) 


Turing instability occurs when E is linearly stable in the absence of diffusion and there 
exists a value of k?, say k?, such that at least one of (37) is reversed. 
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Theorem 6.2 In presence of self-diffusion, Turing instability cannot occur for E, i.e. 
conditions (31) guaranteeing linear stability of E in the absence of diffusion continue 
to guarantee linear stability in presence of self-diffusion too. 


Proof One immediately recovers that if (31) hold, then (37) hold Vk. oO 


6.3 Cross-diffusion driven linear instability 


In order to discover conditions guaranteeing the occurrence of Turing instability, in 
this Section we consider simultaneous linear self and cross-diffusion. Therefore the 
system describing the evolution of a perturbation (U, V, W) to E is given by: 


ðU ~ 
a aU +aj2V +ai3 W + yAU + y2AV + Fi 


aV 2 
ap. = aU +a V + 493W + yn AU + WAV + y23AW + Fo (38) 


aw M 
E = a32 V + a33 W + 32A V + BAW + F3 


whose associated linear system is 


3X yı y2 0 
> = £°X+ D°AX, with DÊ = | yai y y3 (39) 
0 y2 ¥3 


where X and £? are defined in (34)1.2, while det DE = 10123-72332) —V1221V3 > 
0. Let us remark that y;3 and y3; are null because the cross-diffusion between third 
and first species is neglected, being (5) simple food chain. The dispersion relation 
associated to (39) governing the eigenvalues À in terms of the wavenumber k, is 


A3 — TEP? + Fe) — hf (k?) = 0, (40) 


where 


T° (k*) = (L°) — KD) = I} — kyi +v + y), 


IC (KP) = Ky + vivs + v3 — y12V21 — Y23732) + IÈ 


—k [yi (a2 -s)+y2(a11—s)+y3 (a11 +422) — Y12421 
— 721412 — 732423 — V235], 
(k^) = et(D™) + k'[ai1y2y3 + a22.V1 73 — Yı V28 


— Q11Y23 32 + 413 V21 732 + Y12V21S — Yı (V23S + Y32423) 


(41) 


— (yar + y21412)]—k° [yı a228 — a238) — 2411S +73 (411422 —4 12421) 


+ai3(yo1st+y32421) — a11 (V238 + y32a23) + s (Y12421 + y21412)] + Z. 
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Setting 
ClL=V1Y2 — Vi2¥21 + V1¥3 + ¥3¥2 — ¥3223, C3 = I, (42) 
C2=Y1 (d22—S) +2 (411 —S) + 3 (a11 +422) — Y12421 — Y21412 — Y32423 — V235, 
the following theorem holds. 
Theorem 6.3 If (33) holds together with 
c2 > 0, 4c1c3 < 2 (43) 


then Turing instability occurs. 


Proof Conditions (33) imply that E exists and is linearly homogeneously stable. Since 
I 9 < 0> Tez 0, Vk, in order to find conditions guaranteeing the occurrence of 
diffusion-driven instability, one of the following conditions has to be reversed for at 
least ak € Rt: {1£(k?) > 0, hC (k?) < 0, TE(K)IE (k?) — h(k?) < 0}. From 
(41)2 and (42), one has that If (k?) = cik* — cok? + c3, with c3 > 0, since (33) 
implies I > 0. Moreover, since det DE > 0, it easily follows that cı > 0. Therefore, 
a necessary condition to guarantee that I£ (k?) < 0 for some k? € RY, is (43). Let 


. IIE KP) 
us suppose (43); holds. In view of 75 
c 


2 2 ä ä c2 
obtained at (k2) min = TA and is given by (I£ (k?))min = Saa + c3. Hence (43) 


guaranteed that the minimum of ig is negative, i.e. i (k?) assumes some negative 
value. o 


= 2c1k? — c2, the minimum of If (k?) is 


The inequalities (43) define a region where the coexistence equilibrium E is unstable 
and guarantee that Turing patterns occur. 


7 Numerical investigations 


In this Section, we performed numerical simulations to visualize the spatial distribu- 
tions (in a 2D spatial domain) of the solutions to: 


ðu v 

—= 1 A Av, 

ot u| . swt" wR Yi 

ðv u ew 

oe d A A Aw, 4) 
at oe net J+ v + yo, Au + yz Aw 

ðw 


— = ws (1 — 2) + y3 Aw + y32Av, 
ot v 


under the initial-boundary conditions (6)-(7). In particular, we found different sets of 
parameters satisfying conditions (30) and (43) to obtain patterns formation for prey, 
meso and top-predator populations. Based on [32], an original numerical scheme using 
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Table 1 Sets of parameters adopted in the numerical simulations 


Parameter Description Set 1 Set 2 Set 3 
b Inverse of half-saturation constant of preys 0.0499 0.0499 0.0499 
é Top-predators interference on meso-predators 15 0.135 15 

d Death rate of meso-predators 0.9 0.9 0.9 

e Predation rate of top-predators on meso-predators 1.3423 1.3423 1.3423 
n Inverse of half-saturation constant of meso-predators 0.25 0.25 0.25 

q Conversion rate from preys to meso-predators 0.632 0.632 0.632 
s Intrinsic growth rate of top-predators 1 15 15 


the finite difference method for both time and space variables is implemented to explore 
the effect of cross-diffusion on the involved biological species. Table 1 summarizes 
the biological parameters adopted for different simulation sets. In all cases, the final 
simulation time is fixed at 200 days and the required time step interval for discretization 
is At = 10~4. 

With the simulation Set 1 (see Table 1), we observed the almost complete extinction 
of top-predators and a complementary distribution of prey and meso-predator popula- 
tions in the spatial domain. In particular, Fig. 1 shows that the periodical distribution 
of meso-predators is strictly connected to the distribution of preys, which are charac- 
terized by a higher diffusion coefficient (i.e., yj > y2). The occurrence of the Turing 
instability leads to the special case of spotted distributions of meso-predators, which 
are able to further influence the low population density of top-predators. Indeed, the 
top-predators distribution shows higher population densities in the regions where the 
accumulation of meso-predators occurs (see Figs. la and b). 

Successively, we increased the intrinsic top-predator growth rate s and reduced the 
top-predators interference on meso-predators growth rate c in the simulation Set 2, 
see Table 1. By mean of Set 2, we obtained Fig.2. A well-organized distribution of 
preys and meso-predators is clearly obtained in the whole spatial domain. In addition, 
the formation of a structured and spotted distribution of top-predators is also evident. 
However, the only species with significant population density is represented by the 
meso-predators. 

To achieve a more significant scenario from a biological point of view and repro- 
duce a more favourable scenario for prey, meso and top-predator populations, in the 
simulation Set 3 we increased both the intrinsic top-predator growth rate s and the 
top-predators interference on meso-predators growth rate c, see Table 1. 

In Fig. 3, the distribution of prey and meso-predator populations is depicted as a 
self-organized community, but the top-predator distribution is still not clearly periodic 
in space. By increasing the diffusion coefficient y3 and the cross-diffusion coefficient 
y32, Which positively affects the top-predators in terms of final population density, 
we achieved the complementary Turing patterns, extended to all the three species, 
depicted in Fig. 4. The density of top-predators exhibits a notable increase, coupled 
with enhanced mobility across the entire ecological domain, resulting in structured 
spatial distributions around meso-predators. Likewise, the population density of meso- 
predators is higher in areas where top-predators are less abundant, and their spatial 
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Fig.1 Simulation Set 1 after 200 days. The adopted diffusion coefficients are yj = 1, y2 = 0.05, y3 = 1, 
yı2 = 0.003, y21 = 10.1, y23 = 5, y32 = 0.001. a: Spatial domain [0, 20] x [0, 20]; discretization space 
step Ax = Ay = 0.1. b: Spatial domain [0, 5] x [0, 5]; discretization space step Ax = Ay = 0.025 


a 


a 


w 


N 


(b) 


x10 


-3 
20 u 2 v w x10 
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5 0.15 7 
10 25 10 0.1 10 3 
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Fig. 2 Simulation Set 2 after 200 days. The adopted diffusion coefficients are yy = 1, y2 = 0.5, y3 = 
5, yi2 = 9.009, y21 = 52.3, y23 = 0.001, y32 = 8. Spatial domain [0, 20] x [0, 20]. The adopted 
discretization space step is Ax = Ay = 0.1 
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Fig. 3 Simulation Set 3 after 200 days. The adopted diffusion coefficients are y} = 1, y2 = 0.5, y3 = 
3, yı2 = 0.008, y21 = 48.1, y23 = 0.05, y32 = 2.5. Spatial domain [0, 20] x [0, 20]. The adopted 
discretization space step is Ax = Ay = 0.1 
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v w 
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15 7 
1 1 
0 0 0 
0 10 20 0 10 20 


Fig.4 Simulation Set 3 after 200 days. The adopted diffusion coefficients are yj = 1, y2 = 0.5, y3 = 10, 
yi2 = 0.005, y21 = 69.3, y23 = 0.1, y32 = 15. Spatial domain [0, 20] x [0, 20]. The discretization space 
step is Ax = Ay = 0.1 


w 


N 
= 
kej 

N w 


0 10 20 


arrangement is influenced by prey availability. Lastly, prey species predominantly 
exhibit flock-like behavior, as evidenced by localized peaks in population density 
occurring in regions characterized by lower meso-predator densities. 


8 Conclusion 


In this paper, we analysed a three trophic web in the Leslie—Gower predator-prey 
scheme, taking into account a random movement of the populations, hence, both self 
and cross-diffusion are considered. We proved the boundedness of solutions and the 
existence of an absorbing set in the phase space, and we found sufficient conditions 
which guarantee that the governing system does not admit any positive non-constant 
solution. Moreover, via linear instability analysis of the coexistence equilibrium E 
(when it exists), we proved that in presence of self-diffusion, Turing instability cannot 
occur for E, i.e. conditions guaranteeing linear stability of E in the absence of diffusion 
continue to guarantee linear stability in presence of self-diffusion too. Finally, we found 
sufficient conditions guaranteeing the occurrence of Turing instability when both self 
and cross-diffusion are considered and, according to these sufficient conditions, we 
performed numerical simulations in order to visualize the spatial distributions of the 
three populations. 
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